#####
## Replication File, U.S. Data Analyses:
##
## "Toward a Theory and Measure of Racial Benevolence"
##
## Edana Beauvais and Adam M. Enders
##
## Published in Political Behavior, 2026
#####

################################################################################

####
## Load necessary packages
## 
## Analyses conducted in RStudio 1.0.153,
## using R 3.6 (Mac OS)
####

library(readstata13)
library(ggplot2)
library(scales)
library(devtools)
library(dplyr)
library(psych)
library(mirt)
library(MASS)
library(nnet)
library(ggeffects)
library(dplyr)
library(ggplot2)
library(ggthemes)

################################################################################

####
## Open data
####

# Set working directory
setwd("/Users/adamenders/Dropbox/Benevolent Racism - Beauvais & Enders/Data & analysis/Final Files")

# Loading dataset
data <- read.dta13("Clean, U.S. Data.dta")

################################################################################

####
## Figures from Appendix Section S10
####

myvars <- c("police", "adopt", "projects", "segregation", "busing")

# Histograms for appendix
graph_theme <- function (base_size = 18, base_family = "") {
  theme_bw() %+replace%
    theme(aspect.ratio = 1,
          axis.line = element_line(colour = "black"),
          panel.grid.major = element_line(colour = "#d3d3d3"),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(), panel.background = element_blank(),
          axis.title.x = element_text(size=16),
          axis.title.y = element_text(size=16),
          axis.text.x=element_text(colour="black", size = 14),
          axis.text.y=element_text(colour="black", size = 14))
}

data$police_factor <- factor(data$police, levels = c(1, 2, 3, 4, 5),
                             labels = c("Disagree\nStrongly", "Disagree\nSomewhat", "Neither", "Agree\nSomewhat", "Agree\nStrongly"))

tiff("police-US.tiff", width = 7, height = 6, units = "in", res=300, compression = 'lzw')
ggplot(data = data) +
  geom_bar(mapping = aes(x = police_factor, y = (..prop..)*100, group = 1), stat = "count", colour = "black", fill = "grey") +
  xlab("\nTo keep blacks safe, there should be more \n funding for police in urban areas.") +
  ylab("Percentage") +
  graph_theme()
dev.off()


data$adopt_factor <- factor(data$adopt, levels = c(1, 2, 3, 4, 5),
                             labels = c("Disagree\nStrongly", "Disagree\nSomewhat", "Neither", "Agree\nSomewhat", "Agree\nStrongly"))

tiff("adopt-US.tiff", width = 7, height = 6, units = "in", res=300, compression = 'lzw')
ggplot(data = data) +
  geom_bar(mapping = aes(x = adopt_factor, y = (..prop..)*100, group = 1), stat = "count", colour = "black", fill = "grey") +
  xlab("\nMore whites should do their part by \n adopting or fostering black children.") +
  ylab("Percentage") +
  graph_theme()
dev.off()


data$projects_factor <- factor(data$projects, levels = c(1, 2, 3, 4, 5),
                            labels = c("Disagree\nStrongly", "Disagree\nSomewhat", "Neither", "Agree\nSomewhat", "Agree\nStrongly"))

tiff("projects-US.tiff", width = 7, height = 6, units = "in", res=300, compression = 'lzw')
ggplot(data = data) +
  geom_bar(mapping = aes(x = projects_factor, y = (..prop..)*100, group = 1), stat = "count", colour = "black", fill = "grey") +
  xlab("\nHousing projects should be torn down to help \n blacks integrate into safer neighborhoods.") +
  ylab("Percentage") +
  graph_theme()
dev.off()


data$segregation_factor <- factor(data$segregation, levels = c(1, 2, 3, 4, 5),
                               labels = c("Disagree\nStrongly", "Disagree\nSomewhat", "Neither", "Agree\nSomewhat", "Agree\nStrongly"))

tiff("segregation-US.tiff", width = 7, height = 6, units = "in", res=300, compression = 'lzw')
ggplot(data = data) +
  geom_bar(mapping = aes(x = segregation_factor, y = (..prop..)*100, group = 1), stat = "count", colour = "black", fill = "grey") +
  xlab("\nHousing projects make the problem of \n racial segregation worse.") +
  ylab("Percentage") +
  graph_theme()
dev.off()


data$busing_factor <- factor(data$busing, levels = c(1, 2, 3, 4, 5),
                                  labels = c("Disagree\nStrongly", "Disagree\nSomewhat", "Neither", "Agree\nSomewhat", "Agree\nStrongly"))

tiff("busing-US.tiff", width = 7, height = 6, units = "in", res=300, compression = 'lzw')
ggplot(data = data) +
  geom_bar(mapping = aes(x = busing_factor, y = (..prop..)*100, group = 1), stat = "count", colour = "black", fill = "grey") +
  xlab("\nBusing black children to better school districts \n exposes them to values that promote success.") +
  ylab("Percentage") +
  graph_theme()
dev.off()

################################################################################

####
## Factor analysis, reliability, and IRT
####

# Figure S6
check.scree <- fa(data[myvars], fm = "pa", SMC = TRUE, rotate = "none")

tiff("scree-US.tiff", width = 5, height = 5, units = "in", res=300, compression = 'lzw')
xyplot(check.scree$values ~ 1:ncol(data[myvars]),
       aspect = 1,
       type = "b",
       col = "black",
       xlab = "Factor",
       ylab = "Eigenvalue",
       pch = 16
)
dev.off()

# Parallel analysis
fa.parallel(cor(data[myvars]), n.obs = nrow(data[myvars]), fa = "pc")

# Exploratory factor analysis results
fac1 <- fa(cor(data[myvars]), nfactors = 1, fm = "pa", rotate = "none")
fac1

# Checking Cronbach's Alpha
psych::alpha(data[myvars]) 

# IRT for item and test information
grm.mod <- mirt(data[myvars], 1, itemtype = "graded", SE = TRUE)
coef(grm.mod, IRTpars = TRUE)

M2(grm.mod, type = "C2")

plot(grm.mod, type = 'info')

# Figure 2
tiff("info-US.tiff", width = 6, height = 5, units = "in", res=300, compression = 'lzw')
plot(grm.mod, type = 'infotrace', main = "")
dev.off()

################################################################################

####
## Figure 4
####

hist1 <- histogram(~brscale, 
          data = data,
          aspect = 1,
          xlab = "Racial Benevolence",
          ylim = c(-1, 30),
          panel = function(...){
            panel.histogram(...)
            panel.abline(v=mean(data$brscale), lty=2)
          }
          )

hist2<- histogram(~raceresent, 
          data = data,
          aspect = 1,
          xlab = "Racial Resentment",
          ylim = c(-1, 30),
          panel = function(...){
            panel.histogram(...)
            panel.abline(v=mean(data$raceresent), lty=2)
          }
)

tiff("hist-US.tiff", width = 7, height = 3.75, units = "in", res=300, compression = 'lzw')
print(hist1, position = c(0, 0, .5, 1), more = TRUE)
print(hist2, position = c(.5, 0, 1, 1), more = FALSE)
dev.off()

################################################################################

####
## Figure 6
####

# Donation multinomial logit
data <- data %>% 
  mutate(donation = factor(donation, levels = c(7, 8, 10), labels = c("empower", "assimilate", "none"), ordered = TRUE))


mod_donate <- multinom(donation ~ brscale + raceresent + pid + ideo + 
                         college + age + income + female + south, data = data)
summary(mod_donate)

# p-values
z <- summary(mod_donate)$coefficients/summary(mod_donate)$standard.errors
p <- (1 - pnorm(abs(z), 0, 1)) * 2
round(p, 3)

tiff("donateBR-US.tiff", width = 7, height = 3.5, units = "in", res=300, compression = 'lzw')
mod_donate %>% ggeffect(terms = "brscale")%>% #use [all] to get a smooth result
  ggplot(., aes(x, predicted), facet = TRUE) +
  geom_line(aes(color = response.level)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = response.level), alpha = .25) +
  # facet_grid(.~predicted) +
  facet_wrap(. ~ response.level, dir = "h") +
  ylab("Predicted Probability") +
  xlab("Racial Benevolence") +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5),
        panel.grid.minor.x=element_blank(),
        legend.position = "none",
        aspect.ratio = 1,
        axis.line = element_line(size=1, colour = "black"),
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size=12),
        strip.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14, face="bold"),
        axis.title.x = element_text(size = 14, face="bold"))
dev.off()

tiff("donateRR-US.tiff", width = 7, height = 3.5, units = "in", res=300, compression = 'lzw')
mod_donate %>%
  ggeffect(terms = "raceresent") %>% #use [all] to get a smooth result
  ggplot(., aes(x, predicted), facet = TRUE) +
  geom_line(aes(color = response.level)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = response.level), alpha = .25) +
  # facet_grid(.~predicted) +
  facet_wrap(. ~ response.level, dir = "h") +
  ylab("Predicted Probability") +
  xlab("Racial Resentment") +
  #scale_x_discrete(labels=c("empower" = "Educational\nEmpowerment", "assimilate" = "Educational\nImmersion\n(Assimilative Policy)",
   #                         "none" = "None")) +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5),
        panel.grid.minor.x=element_blank(),
        legend.position = "none",
        aspect.ratio = 1,
        axis.line = element_line(size=1, colour = "black"),
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size=12),
        strip.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14, face="bold"),
        axis.title.x = element_text(size = 14, face="bold"))
dev.off()

################################################################################

####
## Figure S4
####

dm <- read.csv("ANES Compare.csv")
dm$var <- reorder(dm$var, dm$est)

tiff("compare-US.tiff", width = 5, height = 5, units = "in", res=300, compression = 'lzw')
ggplot(dm, aes(y = est, x = var, ymin = lb, ymax = ub, color=as.factor(survey))) +
  geom_pointrange(size=.35) +
  coord_flip() +
  theme_bw() +
  scale_x_discrete("Variable") +
  scale_y_continuous("Mean", limits = c(0, 1)) +
  theme(aspect.ratio = 1.65, legend.title=element_blank())
dev.off()

